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Existence of shear horizontal (SH) surface waves in 2D-periodic phononic crystals with an asym- 
metric depth-dependent profile is theoretically reported. Examples of dispersion spectra with band 
gaps for subsonic and supersonic SH surface waves are demonstrated. The link between the effective 
(quasistatic) speeds of the SH bulk and surface waves is established. Calculation and analysis is 
based on the integral form of projector on the subspace of evanescent modes which means no need 
for their explicit finding. This new method can be extended to the vector waves and the 3D case. 
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I. INTRODUCTION 

Emergence of phononic crystals has reinforced the in- 
terest to surface waves in periodic media and heighten 
the need for efficient methods of calculating their dis- 
persion branches. Considerable work has been done for 
solid structures which are periodic along the surface but 
uniform along the depth direction. The latter implies 
pure exponential dependence on depth coordinate and 
thus facilitates the plane-wave expansion (PWE) which 
acts on the surface coordinates only. Applying PWE 
to the wave equation provides a formally infinite alge- 
braic system whose truncation enables explicit finding of 
the evanescent (decreasing with the depth) modes which 
are then used to satisfy another formally infinite sys- 
tem obtained by PWE of the boundary condition on 
the free surface. The above two-step PWE procedure 
was first implemented for Rayleigh waves in a periodic 
structure of layers normal to the surface^ and then ex- 
tended to 2D phononic crystals composed of elastic!^ 
and piezoelectric^' rods normal to the surface. Surface 
waves in such periodic structures uniform along the depth 
direction were also calculated by FDTE03 and waveleP 
methods. 

By contrast, much fewer results for surface waves 
are available in the alternative case of structures which 
are periodic both along the surface and the depth co- 
ordinates, such as rods parallel to the surface. Wc 
have found only three references reporting calculation 
of surface-wave dispersion in depth-dependent phononic 
crystals^ES, both papers using pure numerical means 
(namely, the supercell approximation approach). One 
of the apparent difficulties due to depth dependence 
is a numerically more involved procedure of identify- 
ing of evanescent wave harmonics. Note that the so- 
called extended PWE was suggested as a tool for this 
purpose^HII] but its application to the surface wave 
problem in hand has not been envisaged. 

The present paper pursues two objectives. The first is 
a new method for calculating surface wave branches in 
depth-dependent phononic crystals. The main point and 
advantage of the method is that the dispersion equation 
is expressed in terms of the projector on the subspace of 
evanescent modes and this projector is defined directly 
from the material coefficients expanded in Fourier series 



in surface coordinate(s), without a need to solve for par- 
tial modes and to explicitly sort out the evanescent ones. 
In principle, the proposed method can be used for general 
case of vector waves in arbitrary periodic solid structures, 
but here it is applied to shear horizontal (SH) waves. 
Study of SH surface waves in 2D phononic crystals is 
the second objective of the paper. This is an interest- 
ing problem of its own right. It is well-known that the 
SH surface (localised) waves in ID periodically layered 
half-space can or cannot exist if the layers are parallel or 
orthogonal to the surface, respectively. At the same time, 
we are unaware of results providing explicit evidence of 
SH surface waves in 2D periodic structures. Uncoupling 
of SH modes implies 2D depth-dependent structures, the 
case which defies PWE and was treated by the supercell 



method inSHlil however, the surface waves with SH po- 
larization appear to be beyond the scope of this method. 
As pointed out irP, the difficulty came from the fact that 
the SH surface waves could occur only if the unit cell 
was asymmetric but this made the supercell method " in- 
sufficient or inappropriate". In contrast, asymmetry of 
periodic profile causes no inconvenience for our method. 
By its means we demonstrate examples of subsonic and 
supersonic dispersion branches of SH surface waves in 2D 
phononic crystals. 

The paper is organised as follows. The statement of the 
problem is outlined in Sec. II. The method for calculating 
the surface-wave spectrum is developed in Sec. III. Its 
application is exemplified in Sec. IV. Properties of the 
SH surface waves and the generalization to 3D case are 
discussed in Sec. V. Main findings are summarized in 
Sec. VI. 



II. STATEMENT OF THE PROBLEM 

Consider SH surface waves in a 2D periodic half-space 
{x = (xi,X2) : xi > 0} with a traction-free surface 
X2 = 0. The problem consists of the wave equation 
complemented by the boundary and radiation conditions, 
namely 



'<9i(/z<V) 



■<9 2 (m<9 2 w) = 

0, lim v 

x 2 — yoc 



-pui 2 v, 



(1) 



2 



where di = d/dxi and the shear coefficient p{x) and 
density p(x) are 1-periodic in x\ and x-2. (All subsequent 
results remain explicitly valid for rectangular lattices and 
can be adjusted straightforwardly to the case of oblique 
lattices.) Applying PWE in surface coordinate xx, i.e. 
inserting ID Floquet condition along with the Fourier 
expansion 

u(x)e lklXl with u periodic in xx, 



u(x) 



1 for h 



U,fi,p, 



casts ([T]) in the form 

( (d + ki)n{d + ki)u - d 2 (p,d 2 u) = poj 2 u, 



d 2 u\. 



lim u : 

X2— >oo 



(2) 



(3) 



with 



h(x 2 ) = {h n - m (x 2 )) for h = fi,p, 
d = 2ir(mS n - m ), u(x 2 ) = (u m (x 2 )), n,me1. (4) 

For practical use, we assume all objects in Q to be of 
finite dimension. Equation ^ can be rewritten as 



|Y - Qi?, 

117(0) - (u 0) T , lim rj(x 2 ) = 0, 



(5) 



where ' = d 2 and 

Q 



u 
/xu' 



2 

U) p 



(6) 



The solution to r/' = Qr) with initial data 77(0) is 

ri{x 2 ) = M(x 2 )r]{0) with M(x 2 )= / (I+Qdxa), (7) 

Jo 

where J is the multiplicative integral and I is identity 
matrix. Introduce the monodromy matrix 



Mo = Mo(w,fci)=M(l) 



and let q and w be its eigenvalues and eigenvectors. Tak- 
ing some w as initial data in ([7| defines the Floquet mode 
rj(x 2 ) = M[x2,0]w = q n M.[x 2 — n, 0]w (where n = L^J), 
which is either propagating or increasing or decreasing at 
x 2 — > 00 depending on the absolute value of the eigen- 
value q corresponding to w. In the case of propagating 
modes (\q\ = 1), the spectrum u = w(k) (k = (ki,k 2 ) e 
R 2 ) defined by the equation q(oj, kx) = e lk2 is called the 
Floquet spectrum. 



III. PROJECTOR-BASED METHOD FOR 
CALCULATING SURFACE WAVE SPECTRUM 

1. Projectors. Partition the eigenspace of Mo into 
the folowing subspaces: 



Cd — ( w : M w = qw, \q\ < 1) 
C p = (w : M w = qw, \q\ = 1) 
C x = (w : Mqw = qw, \q\ > 1). 



(9) 



where (...) means a span. Taking ^ with 77(0) from Cd or 
C p or Ci leads to decreasing or propagating or increasing 
solution r](xi), respectively. The projectors P Q on C a 
(a = d,p,i), i.e. 



P„w 



w, we£ Q , 

0, w^£ Q , 



(10) 



can be defined by the formulas 



f p d=2^ / (^I-Mo)- 1 ^, 

|^|=l-0 



2tti ( / 



/ )(zI~M Q )- 1 dz 



Pi = 1 

By definition 



|z|=l+0 1*1=1-0 
1 

27TJ 

|z|=l+0 



(11) 



/ (zI-Mq)- 1 ^. 



Pi 



(12) 



Note that both P p and P; can be expressed in terms of 
P d , see ([381). 



2. Dispersion equation. The surface wave problem 
([5]) is equivalent to any one of the following conditions 



3u ^ 



3u ^ 



u 



e c d 



where 



3u ^ : P d 3U = 0, Pdiu = u , 



Pdl Pd2 
Pd3 Pd4 



(13) 

(14) 
(15) 

(16) 



According to (15), the surface wave spectrum w saw (fci) 



(8) can be defined by the dispersion equation 



D saw (w, fci) = det((P d i - I)* (Pdi - I) + Pd 3 Pd 3 ) = 0, 

(17) 

where * means Hermitian conjugation. 

3. Projection of the Floquet spectrum on the 
plane (lj, k%). Recall that the Floquet spectrum w(k) 
can be determined by the equation 



3w^0: M (o/,fci)w = e lk2 w. 



(18) 



Introduce the multiplicity of the projection of the Floquet 
spectrum on the plane (u>, k\) 



N p (uj, kx) = #{fc 2 G R : w(k) = oj}, 



(19) 



which indicates the number of propagating modes. By 
((9j), ([10]) and ([18]), it can be evaluated as 



N p (u, kx) = dim£ p = Trace P p . 



(20) 



By definition, N p is a piecewise constant function with 
integer even values. Denote the boundaries of areas 
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where N p takes the same value by wtr(fci) and call them 
transonic curves. These curves coincide with the pro- 
jection of local extrema of Floquet branches w(k) on 
the plane (co,ki). The areas of (w, fci)-plane where 
N p ^ 0(<S> P p ^ 0) and A p = 0(<^> P p = 0) will be 
referred to as propagative and non-propagative domains, 
respectively. The whole subsonic range {(u), k\) : uj < 
u>tr,i(ki) = minfc 2 w(k)}, i.e. the part of the plane (co, k\) 
below the minimal frequency of the Floquet spectrum 
w(k), is always non-propagative by definition. Owing to 
existence of spectral gaps in w(k), the non-propagative 
domains may also arise in the supersonic range, i.e. above 

Wtr,l( fc l)- 

4. Refined procedure in the non-propagative 
domains. Consider a non-propagative domain 
Np(w,ki) = 0. According to ( fl5| ) and the identity 
Pd + P; = I (see ([l2])), any root (wjfci) of the equation 

(21) 



D d3 (u,kx) = detP d3 = 



corresponds to a surface- wave solution w saw (fci) (see ( 13 )) 
or a nonphysical solution Wj(fci) which cont ain s of in- 
creasing modes (i.e. (u 0) T e d instead of (13)). Note 
that Dd3(w, fci) is a real function for real arguments, since 
Pd3 is a self-adjoint matrix at N p = 0, see (39). Un- 
like non-negative D saw in (17), the function D d3 gener- 



ally changes sign at its zeroes. Thus seeking the surface 
waves in non-propagative domains, it is convenient to use 
(21 ) alongside (fl7|: the former verifies that a "numerical 



zero" is not a deep but nonzero minimum and the latter 
checks out whether this zero defines a surface wave rather 
than a non-physical wave. Note that surface waves can 



also be identified by using (21 1 along with the condition 
PdiUo 7^ 0, where u is defined by Pd3U = 0. 

In the propagative domains N p (uj,ki) ^ 0, due to 



rankPd = rankP; and P p ^ it follows from (j 1 2|) that 



rank Pd3 < rank Pd < \ rank I = dim Pd3 . Therefore 



(21 ) is an identity in propagative domains and so it can- 
not be used for defining w sa w(&i) in propagative domains. 

5. Calculation of projectors. Equation (11) de- 
fines the projector P d as an integral of the resolvent 
R-zO = ( Z I — Mo) -1 . One way to obtain R 2 o is to com- 
pute M defined by §f§. However, computing M of large 
algebraic dimension 2d x 2c? (which means taking into 
account many members of the Fourier series, see Q) is 
numerically troublesome because some of its components 
grow exponentially as d increases. On the other hand, 
components of R z o can grow in general only as fast as 
linearly in d. Therefore it is numerically adavantageous 
to calculate R z o directly rather than via Mo. 

Denote 



R«(x 2 ) = (al - M{x 2 )Y 



(22) 



where a is some fixed constant which does not belong to 
the spectrum of M. As a solution to (5]), M satisfies the 
linear differential equation with initial data 



M'(x 2 ) = Q(x 2 )M(x 2 ) 
M(0) = I, 



(23) 



from which it follows that R„ satisfies the Ricatti equa- 
tion 



R-'aO*^) = H a (x 2 )Q(x 2 )(a'R a (x2) 
R Q (0) = (a- l)- 1 !. 



(24) 



Numerical integration of (24) provides the value R„o 



R Q (1) for a given a. Once it is found, the identity 



R~ 



R a0 (I+ {z - a)R tt o) 



(25) 



yields R z0 with varying z as required in ( |11| . Thus com- 
bining (11) and (25) defines the projector Pd in the fol- 



lowing specialized form 



RaO 



re llp {l + {re iv - a)R a0 ) -1 ^ (26) 



with r < 1 sufficiently close to 1. Two other projectors 



Pi and P p are expressed via Pd, see (38). 

In numerical implementation, we used the fourth or- 
der Runge-Kutta method for calculating R Q o from (24) 
where a 6 C has being chosen randomly anew for each 
next calculation point (lj, k\). It is recommended to 
avoid taking a near the real axis or the unit circle (we 
used a £ [—6, —3] x [3i, 6i]). The Chebyshev method was 
employed for evaluating ( [26] ) where we took r = 0.99. 
Dealing with Eqs. (17) and (21) it is convenient to use 



'normalised' functions 

£> saw (u;,fci) = A min (w,fci), 
D M (u,ki) = det(P d r 3 1 (0,fc 1 )Pd 3 (w,fci)), (27) 

where A m j n is the minimal eigenvalue of the self-adjoint 
non-negative operator whose determinant is taken in 
{17|). 



IV. EXAMPLES 

1. Stiff cylinders in a soft matrix. Assume an 
epoxy matrix with a periodic structure of steel cylindri- 
cal bars parallel to the free surface x 2 = 0. The material 
constants are p — 1.14 g/cm 3 , fj, = 1.48 GPa for epoxy 
and p = 7.8 g/cm , p = 80 GPa for steel. Let the steel 
bars form a rectangular lattice with the horizontal and 
vertical periods ai : a 2 — 1 : 2, and the radius of bars be 
0.45ai (in the following, a% = 1 and k\a\ = ki, Loa\ = u 
in units mm/|is). For obtaining SH surface waves, it is 
necessary that the unit cell is asymmetric about the hor- 
izontal midplane (see Sec. V.l). We shall consider two 
interrelated types of structures which are reciprocal to 
each other in the sense that their unit cells pass into one 
another by means of reflection about the horizontal mid- 
plane, see Fig. 1. Note that the unit cell of the structure 
(a) is 'faster' on the top than on the bottom, hence vice 
versa for the structure (b). Figure 2 demonstrates that 
the two reciprocal structures are characterized at fixed 
(lj, ki) by the same number N p of propagating modes, 



4 



whereas the frequency w saw of surface- wave solutions for 
the configuration (b) is the frequency U{ of non-physical 
solutions for the configuration (a). A general proof for 
this feature is given in Sec. V.2. 

The surface- wave dispersion branches w saw (fci) for the 
structure (b) are shown in Fig. 3. For this calculation, 
we used d = 17 terms of the Fourier series ©. Nu- 
merical data clearly show that there exists a subsonic 
branch w saWj i below the first transonic curve w tr ,i, though 
their relative difference is only of the order of 1%. Sub- 
sonic and supersonic surface waves occurring in the non- 
propagative domains (N v = 0) are defined by common 
zeros of Eqs. (17) and (21 1. In propagative domains 
(N p ^ 0), a surface wave is defined only by (17) and 



AD s&w ,D d3 ,N p 



is therefore associated, in the numerical context, with a 
minimum which tends to zero with growing number d of 
terms of Fourier series. The surface- wave branches shown 
within the upper propagat ive domain in Fig. 3 yield the 
value of -D s;i , 



Amin (see (27)) of about 10 



stress free 



stress free 




FIG. 1: Reciprocal configurations of a periodic structure of 
steel cylinders in epoxy matrix. 



2. Perturbation of ID-periodic structure. In 

order to illuminate formation of surface-wave spectrum, 
let us think of a vertically periodic stack of equidistant 
epoxy and iron layers (Fig. 4a) and then assume that 
each epoxy layer contains narrow lead plates embedded 
periodically along the horizontal direction (Fig. 4b). The 
material constants arc p = 116 g/cm 3 , p = 7.88 GPa for 
Fe and p = 14.9 g/cm 3 , p = 11.6 GPa for Pb. Taking the 
ratio of the horizontal to vertical periods a\ : a 2 = 1 : 3 
for both cases and applying the calculation procedure of 
Sec. Ill leads to the spectra shown in Figs. 5a and 5b, 
respectively. It is clearly seen that introducing horizontal 
periodicity creates band gaps at the edge and inside of 
the Brillouin zone. 




UJi O UJt 




FIG. 2: (color online) The function N p (uj) (see pM) and 
the normalized functions -D S aw(<^) and Dd3(w) (see ( 27 1 ) cal- 
culated at fixed fci = 4-7r/5 for the reciprocal structures (a) 
and (b) of Fig. 1, respectively. The curves Dd3 and N p do 
not change on passing from one configuration to another, see 



(31 1 . Zeros of D^3 at jV p = are surface- wave solutions u) B: 



when they coincide with zeros of -D S aw, 
non- physical solutions oJi. 



otherwise they are 



V. DISCUSSION 

1. Reciprocity property for SH surface waves. 

Given periodic functions /i(x) and p(x) of x = (x%, £2) € 
M 2 , let us mentally cut the space into two halves by the 
plane X2 = and turn the half-space X2 < upside 
down. Thus we obtain two models of a half-space x-i > 
0: one with the profile /i(x), p(x) and another with the 
reciprocal profile /I(x) = p(xi, —x 2 ), p(x) = p{x\, —x 2 ). 
Consider the relation between the properties for a direct 
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k Co>, mm/ fis 




k CJ, mm//j,s 



FIG. 3: Surface- wave dispersion branches oj sa w(fci) and tran- 
sonic curves uj t i(ki) for the structure displayed in Fig. lb. 
Blank and hatched areas are the non-propagative and prop- 
agative domains, respectively. 



stress free 



stress free 




(a) 



(b) 



FIG. 4: Periodic structure of epoxy and iron layers (a) and a 
similar structure but with lead baffles in epoxy (b). Lead to 
epoxy volume fraction is 1 : 9. 




(a) 



. UJ : mm//j,s 




(b) 



FIG. 5: Surface-wave dispersion spectra aj saw (fci) (a) and (b) 
for the structures displayed in Fig. 4a and 4b, respectively. 
Notations are the same as in Fig. 3. 



and reciprocal profiles (the objects constructed from /i, 
p will be labeled by a tilde) . 

By Q and the definition of multiplicative integral, 



N- 



1 

r 



Mo X = Jim n (I- 



N-yoo j =N 



Jim J] (I-^Q(^) ) -SM S (28) 
I 



with S = S 1 = 







where we used that SQS = — Q since Q has zero diagonal 
blocks, see According to (28 1, the eigen-subspaces 
and projectors for direct and reciprocal profiles are re- 



lated as follows: 



C d = SA, Pd = SPjS, 
A = S£d, P; = SP^S, 



(29) 



Hp S£p, 



Pp = SP P S. 



Using a similar to (16) notation for the blocks of Pj, 



introduce the equation 

A(w,*l) = dBt((Pii-I)*(P U 



0, (30) 



whose solutions describe non-physical waves which 

satisfy the stress-free condition but consist of increasing 
modes. From p9l, with reference to (Il7|), pfl and pot, 



D 



d3 



D d3 , N p 



N p . (31) 



G 



It is thus proved that (i) the surface wave solution 
Wsaw(^i) for a direct profile p, p is at the same time a 
"non-physical" solution uj\(k\) for a reciprocal profile p, 
p and vice versa, and that (ii) the number of propagat- 
ing modes for direct and reciprocal profile is the same. 
Also we obtain from ( 29 1 that (iii) there is no SH surface 
waves if the profile p, p is symmetric in the depth direc- 
tion, i.e. if p = fx, p = p. Indeed, assume the opposite: 
there exists a surface wave in the case of a symmetric 
profile. Then using (13 1, (29 1 and p. p = p, p yields 



= S 



e SCd — Ci — £i, 



(32) 



which is in contradiction with Cd (1 C\ = {0}. 

2. Effective (quasistatic) speed of the funda- 
mental SH surface wave. The surface-wave spectrum 
w saw (fci) may or may not contain the so-called funda- 
mental branch which starts at zero uj and fci. Suppose 
that it does, i.e. that the first (lowest) branch w saw l (fci) 
is a fundamental branch. Denote 



lim 

fci— s-0 



c(/c) = lim 



^saw.l(fcl) 

wi(k) 



|kP"o Ikl \ K |k 



(33) 



Ctr 



= lim 



Wtr,l(fel 



fci->o ki 



where c saw is the effective speed of the fundamental SH 
surface wave, c (n) is the effective speed of bulk Floquet 
modes (wi(k) is the lowest sheet of the Floquet spec- 
trum), and Ct r is the onset slope of the first transonic 
curve wtr,i(fci)- It can be proved that 



Ctr- 



(34) 



On the other hand, from ( 33 1 and the definition 
w tr,i(fci) = mm fc 2 w i(k) we deduce that c 2 r — 
min^2 k^ 2 c 2 (/c). Recalling that a squared c(/c) is a 
quadratic form of k, it follows that 



1 



det (Cij) , where c 2 (k) — KiCijKj (i,j = 1, 2). 

(35) 

Explicit expressions for the components of the matrix 
(Cij) can be found in Sec. IIC oP^. Note that this matrix 
is diagonal for the structures displayed in Fig. 1 and non- 
diagonal for the structure displayed in Fig. 4b. 

3. Algebraic symmetries. Recall that ^ and Q 
imply the standard identities 



Q 



-T _1 QT, M 1 = T _1 M*T 



where 



I 

-I 



(36) 



(37) 



Using (36) with (11) gives 

Pa = ^ 



J (zl-M*)- 1 ^* 
l«l=i-o 

— / {zI-M^dz 
m J 

|*|=i-o 

".— / (zI-TMoT^dzT* 
2m J 

l*l=i-o 

?i / (zI-M^ 1 )" 1 ^* =T" 1 P i T. 
2m J 



|*|=i-o 



Thus, with reference to (12), knowing the projector Pd 



yields the two other projectors as follows 
p. = x -i p * T 

P p = I P d T : P*T (=T 



l p;t) 



(38) 



In the non-propagative domains N p — 0(<=> P p = 0), Eq . 
Q reduces to P d + T~ 1 P d T = I and so the blocks ( [i~6[ ) 
ofPd satisfy 



P d2,d3 — Pd2,d3, 



dl 



d4 



(39) 



By ( |36| , the ^-component of energy flux averaged over 
x\ and time, V{x%) — r/' 'IV/. is non-zero or zero if, 
respectively, r\(x2) — M(x2)w is a propagating or non- 
propagating Floquet mode, as defined in Sec. II. Equal 
number of positive and negative values of T 3 , i.e. of for- 
ward and backward propagating modes, is due to zero 
signature of T. 

In conclusion, note that the definition of r) in (JeJ) differs 
from another conventional form which involves an addi- 
tional factor V and hence leads to the similar identities 
for Q and M but with T having both off-diagonal blocks 
equal to I (unlike p7|)). 

4. Spectral decomposition in the non- 
propagative domains. The integral definition (11) 
of the projectors underlying the present method on the 
whole and the above identities in particular is indepen- 
dent of accidental degeneracies of eigenvalues q of Mo. 
Let us assume a generic case where all q are distinct and 
hence Mo possesses a full set of eigenvectors w. Then 
besides (11) projectors satisfy the spectral decomposi- 
tion. Consider a non-propagative domain (all \q\ =/= 1). 
Apply the normalization w^Tw; = 1 where Wj and W; 
correspond to |g| < 1 and |g| > 1, i.e. to decreasing 
and increasing Floquet modes. Define the matrix W 
whose columns are sets of normalized Wd = (a d b d ) T 
and Wj = (a; b;) T , namely, 



W 



A d A; 
Bd Bj 



(40) 



where blocks Ad,i and Bd,i consist of ad.i and bd,i- By 
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(36), WT*W*T = I and hence 



(A d 








\B d 


o 




-f) - ( 



A d B* 
B d B* 



—Ad A* 
B d A* 



P; 











G 




V- B 5 





-AB d 

BiB* 



AiA d 
BiA d 



,(41) 
.(42) 



Note that the structure of the left off-diagonal block P d3 
corroborates with the interpretation of Eq. (21 ) given in 
Sec. III.4. 

5. General 3D case and the depth-independent 
case. The projector-based method described in Sec. Ill 
can be extended to vector surface waves (Rayleigh waves) 
and to 3D periodic materials by means of a single replace- 
ment of the SH matrix Q by its general form 



Q 



c- 1 



-pu? -A 2 - A\C~ 
C = (cfui), Ai = (caki(di + iki)), 
A 2 = ((dj + ikj)cij k i(di + ikt)), j, 1 = 1,2, 



(43) 



where Cyjy(x) are the stiffness coefficients. For brevity, 



(43) is written in the x-space and is subject to Fourier 
expansion in the surface coordinate (s). Note that (43) 
with fc 1)2 = reduces to (5.28) oP 1 . 

If Cijki are uniform in the depth direction X3, then 
the definition of projectors P can be simplified due to 



replacing Mq in (11) by 



(44) 



L = (Q + oI)(Q-a*I)" 



with some random a £ C and Re a > 0. Using (|44j) fa- 
cilitates the calculations significantly because it removes 
a need to solve the Riccati equation ( 24 ) for finding the 
projectors. 



SH surface waves in 2D semi-infinite phononic crystals 
periodic along the surface and depth coordinates x\ and 
x 2 . For instance, these may be periodic structures of ho- 
mogeneous bars parallel to the free surface x 2 = 0. The 
necessary condition for the existence of SH surface waves 
is vertical asymmetry (in x 2 ) of the unit-cell profile. By 
analogy with the ID-periodic case of layers parallel to the 
surface (see^j), the subsonic SH surface waves are more 
probable if the unit cell is 'slower' on the top than on the 
bottom. The onset slope of this branch, i.e. the effective 
(quasistatic) speed of SH surface waves, is the slope of 
projection on the plane (u, k\) of the lower bound of Flo- 
quet spectrum w(k) of the infinite 2D-pcriodic medium. 
Enhancing the unit-cell horizontal asymmetry (in x%) in- 
creases the band gap which separates the subsonic branch 
and the next surface-wave branch at the edge of Bril- 
louin zone. Turning the profile upside down replaces the 
SH surface wave solutions by the 'non-physical' ones (in- 
creasing into the depth). 

The method proposed herein for calculating the sur- 
face wave spectra is based on the dispersion equation ex- 
pressed through the projector on the subspace of evanes- 
cent modes. It is defined as an integral of the resol- 
vent of the monodromy (transfer) matrix, and the in- 
tegrand function satisfies Ricatti equation whose coeffi- 
cients are members of ID Fourier series of material prop- 
erties. Thus there is neither a need to identify partial 
modal solutions of the wave equation, nor to calculate the 
monodromy matrix which is prone to numerical instabil- 
ity. Knowing the projector for the evanescent modes also 
yields the number N p of the propagating modes at any 
given ui and k\ . Root finding of the dispersion equation 
can be refined in the non-propagating domains N p = 0. 
The methods can be generalized for the vector waves and 
the 3D case. 



VI. CONCLUSION 

We have studied the frequency versus horizontal 
wavenumber spectrum co(ki) of subsonic and supersonic 
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